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Abstract 



We compare Kohn-Sham results (density, cohesive energy, size and effect of charg- 
ing) of the Spherical Averaged Pseudopotential Model with the Stabilized Jellium 
Model for clusters of sodium and aluminum with less than 20 atoms. We find that 
the Stabilized Jellium Model, although conceptually and practically more simple, gives 
better results for the cohesive energy and the elastic stiffness. We use the Local Den- 
sity Approximation as well as the Generalized Gradient Approximation to the exchange 
and correlation energies. 
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1 Introduction 



The energy functional for an electronic system under the influence of a positive charge 
distribution (ions) is 

E[n, n + ] = T[n) + E xc [n] + \ [ d 3 r I d\' n ^ n ^ (i) 
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where n{r) is the density of the valence electrons, n + (f) is the ionic density, T is the kinetic 
energy, E xc is the exchange and correlation energy, and V + (r) is the potential created by the 
ions. Atomic units are used in this paper. 

The Jellium Model is a simple model to describe simple metals. To study a spherical 
cluster in this model, the ions are replaced by a continuous positive background which is 
constant inside a sphere of radius Rj e i = r s iV 1 / 3 , where N is the number of valence electrons, 
and zero outside: 

n+ (r)=n9(R je i-r) , (2) 
with n = 3/(47rr 3 ). For this model the external potential is simply 
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(r > R jel ). 



The Stabilized Jellium Model (SJM) represents an improvement over the Jellium 
Model. For the bulk matter, it yields more realistic binding energies and bulk moduli [I| . For 
semi-infinite systems, it provides more correct surface energies and work functions For 
clusters, it gives very reasonable cohesive energies, dissociation energies, ionization energies, 
etc. P, P| . Its energy functional may be written as 



E SJ [n,n + ] = Ej[n,n + ] + (5v) / d 3 r9(r - R jel ) [n(r) - n + {r)} + e / d 3 rn + (r) , (4) 



where Ej[n,n + ] is the jellium energy functional (Eq. (1) with the external potential given 
by Eq. d)), and 
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is the average, in the Wigner-Seitz cell (taken to be a sphere of radius tq = r^z/ 1 / 3 , with v 
the valence), of the difference between the pseudopotential 

( \ / (r <r c ) , , 

w[r) = < , ) ( (6) 
I — vjr [r > r c ) 

(Ashcroft empty core potential 0, with core radius r c ) and the jellium potential. 

In the last term of Eq. (|2]), e is the spurious background repulsion inside each Wigner- 
Seitz cell, which is equal to the sum of the Madelung energy 

1 r° , . 2-(~ v 1 \ 9 

cm = ~ drArcr n[ h -v jei = - — — (7) 
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and the average of the pseudopotential repulsion JTJ : 

- [ rc dr4nr 2 n- = - r 4 . (8) 

Physically, Eq. means that we start with jellium and localize uniformly the ions but 
their interaction with the valence electrons (accounted for by the pseudopotential) is only 
taken into account perturbatively and in a spherical averaged way. 

With the energy functional @ we may extract results for clusters, via the Kohn-Sham 
method in the Local Density Approximation (LDA) for exchange and correlation, going from 
the single stabilized jellium atom all the way up to the bulk solid. We use the Perdew-Wang 
interpolative formula for the correlation energy. We consider that the jellium background 
is either fixed with the bulk density value or relaxed in order to give the minimal energy 
for a given number of particles (using a fixed pseudopotential, e.g., transferred from the 
bulk). In the latter case, which is in principle more realistic, we speak of self-compression 
of the neutral cluster, a phenomenon which may be classically explained by the effect of the 
surface tension on the density [0, |, [|. We take here, when not otherwise indicated, the 
SJM including self-compression. The energetics of the cluster ground states is similar to 
that of the Jellium Model, with the same shell structure, but shifted down, i.e., the valence 
electrons are more bound. 

Let E = E(N, r s , r c , v) be the binding energy of a neutral spherical cluster with iV 
valence electrons. The equilibrium ionic density, r*, of a cluster is obtained imposing the 
following equilibrium condition: 



d_ 



E(N,r s ,r c , v) 



0, (9) 



N 

where the derivative is taken at constant N and r c = r c (rj, v), where rf is the bulk density 
parameter. The core radius is fixed by the bulk stability condition, which is Eq. (9) in the 
limit iV —* oo. We have self-compression if r* < rf and self-expansion if r* > rf. 
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The elastic stiffness of the cluster is defined by 

E(N,r*„r e ,u) 



In the limit N — > oo, this quantity goes over to the bulk modulus. Stability demands that 
B > 0. 

On the other hand, the Spherical Averaged Pseudopotential Model (SAPS) takes into 
account the geometrical arrangement of the ions located at positions Rj P~G| , |TT| . This is 
optimized by minimizing the total energy of the cluster. 

Only valence electrons are explicitly treated, as in the SJM, being the ion cores replaced 
by atomic pseudopotentials w(\f — Rj\). The simplification introduced by the SAPS model 
for clusters of simple metals consists of replacing the external potential felt by the valence 
electrons, V+(f) = J2j at w(\r — Rj\), by its spherical average, (V+) (r), around the cluster 
center. In this way, the energy functional in the SAPS model reads as 



B(N,rl 



-V 



d 2 E , 



d 2 



\N 



12nr*N dr 2 s 



EsAPsin, Ri} = T[n] + E xc [n] + J / d 3 r [ d 3 r' n ^<^ (n) 

2 J J \f — r'\ 

+ 47rr 2 / 'dr{V + )n(r) + lj:-J%- , 
J z ijLj \Ri — Rj\ 

where denotes the charge of the ions and N at = N/u is the number of atoms. The last 
term in Eq. (11) represents the repulsion between the ion cores, which is approximated 
by the interaction between point charges. As in the SJM, we have used the LDA with the 
correlation energy functional of Perdew-Wang. 

In order to determine the cluster energy, we start from a randomly initial atomic con- 
figuration Rj, and solve the Kohn-Sham equations, evaluating the valence electron density 
and the total energy given by Eq. (11). To obtain the geometry that minimizes the total 
energy, we calculate the energy variation (SEj) a , corresponding to the displacement of the 
ion at position Rj in each Cartesian direction (a = x,y,z), by a small quantity Sd, keeping 
the electron density and the other ions coordinates frozen. Then we choose the maximal 
energy variation 5E max , and move each ion coordinate by the amount [(8Ej) a /8E max \ 8d, ob- 
taining a new geometry, Rj, for which we calculate the corresponding electron density. With 
this new geometry we repeat the cycle again, being the process iterated until convergence 
is achieved and an equilibrium geometry is found. This generally leads to a local minimum. 
To find the global minimum we have repeated the calculations for a large number (typically 
20) of random initial configurations. 
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To obtain the elastic stiffness in the SAPS model, we use an analogue of Eq. ([TQf) with 
the approximation 

d 2 E E(Rj - Sr) + E{Rj + Sr) - 2E{Rj) 



dr 2 \Sr\ 



(12) 



The energy of the cluster is calculated with the atoms at the equilibrium positions Rj and 
displaced radially in and out by a small amount, Sr, from Rj (for simplicity we have omitted 
the energy functional dependence on n). 

An extension of the SAPS model which exploit all ionic degrees of freedom in three di- 
mensions but restricts the electronic many-body problem to axial symmetry has been recently 
introduced by Montag and Reinhard [|T^, O . This model is known as CAPS (Cylindrical- 
Averaged Pseudopotential Scheme), since a cylindrical average is taken instead of a spherical 
average. 

In this work we compare the SAPS model with the SJM because the two schemes use a 
spherical average of the atomic pseudopotential, and one could expect the two descriptions 
of some clusters properties to be similar. In particular, one could expect that the phenomena 
of self-compression of neutral clusters with respect to the bulk density and self-compression 
or self-expansion of charged clusters, which have been identified in the SJM || [|, are also 
described in the SAPS. With this goal, we examine the two models for small clusters of 
sodium and aluminum. In the Appendix we compare the LDA cohesive energies with the 
Generalized Gradient Approximation (GGA) ones. 



2 Results 

Fig. 1 shows the SAPS valence electron density for 8-atoms clusters of sodium (r s = 3.93, 
v = 1) and aluminum (r s = 2.07, v = 3), in comparison with the SJM one, if the jellium 
edge is allowed to move. In the SAPS, we have used the recently proposed evanescent core 
pseudopotential model, which has the advantages of being local and analytical and which 
describes generally well the main physical properties of bulk simple metals |14|, [15], |16| . In the 
SJM, we have used the simpler and more common Ashcroft pseudopotential with the core 
radius fixed by the condition of bulk stability. We could as well have used any other local 
pseudopotential with a single free parameter, arriving at the same results. We conclude that 
the two densities are somewhat different, especially in the interior for aluminum (the SAPS 
density is much lower there). The surface diffusivity is bigger for SAPS than for SJM (this 
leads to a larger redshift of the surface plasmon resonance with respect to the Mie peak |T7| ). 
However, the radial probability density r 2 n(r), represented in the inset, is similar in the two 
models. We have performed the same comparison for positively ionized systems (Nag + and 
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Alg 2+ ), arriving at the same conclusions. 

Let us consider the energy of the cluster in the two models. The cohesive energy is the 
difference between the energy of the free atom and the energy of the cluster per atom, 

E coh (N a t) = E{v) - . (13) 

For the bulk, Eq. flTB"| ) leads to the heat of sublimation: E co h(N at = oo) = E(y) — a v u, where 
a v is the bulk energy per valence electron. Fig. 2 represents the cohesive energy for small 
clusters of sodium and aluminum using the SAPS model, and the two versions of the SJM: 
one with the self-compression effect and another with the background density frozen at the 
bulk value, rf. In Fig. 2 we have also plotted the liquid drop energy in the SJM 

Et D h M {N at ) = a s v^ (1 - iV7 t 1/3 ) + a c *//3 (1 - N^) , (14) 

with a s and a c the surface and curvature energy coefficients. In the case of self-compression, 
the curvature energy coefficient is replaced by 

~a c = a c -l^f. (15) 

For sodium we have a c = 0.18 eV, and for aluminum a c = —0.10 eV, which contrast with 
the non-compression values, respectively 0.26 eV and 0.65 eV. The surface energy (a s = 0.57 
eV for sodium and 0.86 eV for aluminum) is not affected by compression. 

The cohesive energy in the SAPS is much lower than any of the SJM results, although a 
similar pattern (with maxima at shell closures) is apparent. For sodium, the self-compressed 
SJM results are in very good agreement with LDA Car-Parrinello calculations |18| (and also 
close to the experimental values For aluminum, the frozen density SJM results are 

in good agreement with Car-Parrinello calculations [^0] (and experiment pi]), while self- 



compressed SJM and SAPS yield negative values for very small clusters. This fact should 
be seen as a drawback of the spherical approximation used in both models. The spherical 
restriction, valid for the atom, is indeed an artificial constraint for the bulk solid. This makes 
the SAPS valence-electron energies too high and the cohesive energies too low. In the SJM, 
the errors made for the atom and the solid seem to be similar. 

The SAPS result for the cohesive energy depends somewhat on the selected pseudopo- 
tential, but replacing the evanescent core pseudopotential by a different one such as Manni- 
nen's p2| , the SAPS results remain negative for small clusters. The error made in using the 



pseudopotential for the atom may be estimated replacing the binding energy of the pseudo- 
atom by the sum of the first v ionization potentials in an all-electron calculation using the 
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same exchange and correlation energies. That correction, using the evanescent core pseu- 
dopotential, yields a shift upwards in the cohesive energy of 0.2 eV for sodium, and 0.0 eV 
for aluminum. 

In Appendix, we compare the LDA with the GGA results for the cohesive energy, 
concluding that the density gradient corrections reduce this quantity. 

Besides the electronic shell structure, leading to magic numbers, which appears in both 
the SJM and the SAPS model, there is in the latter an ionic structure which also reflects 
itself in oscillations of the total energy. In order to disentangle the electronic and geometric 
shell effects, we have evaluated the kinetic energy contribution by means of an Extended 



Thomas-Fermi energy functional with fourth-order gradient corrections (TFDGW4) |23|, ^4 



using, however, the self-consistent Kohn-Sham density of each model (Fig. 3). In this way 
the electronic shell-structure due to the quantum kinetic energy operator of the Kohn-Sham 
equations is partially erased. In both cases, it is clear that with the new functional we 
would obtain self-consistently a different ionic background and ionic structure, but we have 



just evaluated a semi-classical kinetic energy with the quantal density as in Ref. |24j|. The 
resulting binding energies of sodium neutral clusters are represented in Fig. 3. In the SJM, 
we obtain a smooth curve, which may be fitted by a liquid drop formula; the surface and 
curvature coefficients are a s = 0.52 eV and a c = 0.22 eV, which agree very well with results 
obtained from the planar surface problem, respectively 0.58 eV and 0.26 eV |2|. But, in the 
SAPS, we are left with some wiggles, which are mainly due to the reorganization of ionic 
shells. An increase of the energy arises when a new shell is added. For instance, this happens 
when going from N at = 4 to N at = 5 or from N at = 13 to N at = 14. Note that, except for 
the N,a — 1 case, the SJM results is always above the SAPS results. 

For the sake of comparing the cluster size evolution in the two models, we have plotted, 
in Fig. 4, the radius of the outer ionic shell R out in SAPS divided by R(N) = Rj e i(N) — d/2, 
where d is the distance between parallel closed packed planes in the bulk solid, against 
the number of atoms. In fact, the planar surface is the limit of the curved surface of a 
big cluster, and it is known from surface physics that the first lattice plane is located at a 
distance x = —d/2 from the jellium edge located at x — 0. The quantity d = 1.436 Tq is the 
distance between the (110) planes in the bcc structure (these are the planes with the biggest 
separation) and d = 1.477 Tq is the distance between the (111) planes of the fee structure. 
We see that R ou t/R approaches 1 from above, although it becomes increasingly difficult to 
apply the SAPS to very big clusters and to attain the planar surface limit. This dilatation 
effect of the outer ion shell with respect to the surface limit is particularly clear for very 
small clusters. On the other hand, we observe contraction in the SJM: in Fig. 4, r*/rf goes 
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from below to 1 when N at increases. In both models, the approach to the asymptotic limit 
is more rapid for sodium than for aluminum. 

Notwithstanding the expansion of the outer ionic shell, the distances between next 
neighbour ions inside each shell are smaller than in the bulk environment differing from ionic 
shell to ionic shell. This inhomogeneous contraction observed in SAPS has been reported 



before |25|, ^6[. Such an effect is outside the scope of the SJM, since the jellium background 
must have a constant density [p]]. 

To analyse the behaviour of charged systems we may keep the size fixed (e.g., Nag and 
Alg) and increase the charge all the way up to Coulomb explosion of the volume (Fig. 5). 
In the first case, the radius of the jellium sphere is increasing, while in the SAPS the same 
happens with the radii of the ionic shells. We see a strong similarity of the phenomenon of 
volume explosion described by the SJM and the SAPS model. The maximal positive charge 
that Na§ can hold is q = 3, whereas for Al 8 is q = 13 for the SJM and q — 14 for the 
SAPS (the charged systems of Fig. 1 are therefore close to their stability limit). For Alg", 
we observe in both models a small self-compression instead of self-expansion with respect to 
the neutral system. This is an interesting effect since it goes against the expectations. 

The type of instability we are describing is different from usual fission since, in the 
latter, volume is fixed while shape is deformed and here volume changes keeping the shape. 
The fissibility, defined as x = Ec/2E$, where Eq is the Coulomb and E$ is the surface 
energy of a spherical cluster, controls the fission probability. For x = 1 the fission barrier 
vanishes and the cluster decays spontaneously. A different parameter should be considered 
for the explosion of the background. Anyway, within the SJM, we get x = 3.4 for Nag and 
x = 26.6 for Al 8 as critical fissibilities for volume explosion, which indicates that excited 
clusters prefer to lower their energy by deforming and breaking into two or more pieces than 
by exploding the volume. 

Finally, in Fig. 6 we compare the elastic stiffness of charged aluminum clusters. The two 
values agree when the charge is big, but differ for small charges: the elastic stiffness of neutral 
aluminum clusters in the SAPS model is much larger than in the SJM. We remind that the 
aluminum bulk modulus given by the SJM is a factor of 2 bigger than the experimental 
value, as indicated in the picture. 



3 Conclusions 

We have examined comparatively some physical properties of small clusters of sodium and 
aluminum as given by the SJM and the SAPS model. 
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The radial electronic probabilities are similar for sodium and somewhat different for 
aluminum. Although the SJM is conceptually and practically more simple, its results for 
the cohesive energy are in better agreement with more involved theoretical results. The 
restriction to the spherical symmetry seems to be a serious drawback of the SAPS model. 
The outer ionic radius in the SAPS and the jellium edge in the SJM for neutral clusters 
approach the asymptotic limit in different ways. However, the results of the two models 
compare very well when we normalize the ionic radius of a charged cluster to that of the 
corresponding neutral cluster, and keep the number of atoms fixed when increasing the 
charge. The elastic stiffness given by the SJM, although too high in comparison with data 
for high density metals such as aluminum, is better than the SAPS output. This indicates 
that the spherical restriction to the ionic structure, besides giving too low cohesive energy 
for high density metal clusters, also makes them too hard. Therefore, the corresponding 
monopole compression mode in the SAPS has an energy which is expected to be too high in 
comparison with experiment [^7j . 

In conclusion, the results obtained with the SAPS for the energetics and compression 
properties of small clusters should be taken with some care. Simpler models, such as SJM, 
may surprisingly perform better. On the other hand, the SJM is useless for heterogeneous 
clusters, while SAPS is able to describe them. The two types of models are useful for 
understanding the main trends of small and large clusters, since ab initio theoretical models 
are computationally more elaborated and, at present, almost impossible to apply to clusters 
with more than 20 atoms. Since the SJM and the SAPS model assume spherical symmetry, 
they are particularly suitable for systems with spherical geometries (for instance, N at = 8 
and 20 for sodium). They are not adequate for systems with N at < 5, which are known, 
from first principles calculations, to be planar. 

The assumptions of stabilized jellium are not restricted to a compact spherical shape 
making it more versatile than the SAPS model to describe some properties of metallic clus- 
ters. A simple modification of the SJM which should bring it to a better agreement with 



SAPS is the hollow SJM [^, |29[. Indeed, for a small system, we may open a hole in the 
middle and vary the size of the system so that the energy is minimal for a given number of 
particles. Deformed jellium [B(J and CAPS |12], [13|] represent improvements of respectively 
spherical SJM and SAPS which correct the spherical approximation. 
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Appendix: The Generalized Gradient Approximation 
for small clusters 



The Local Density Approximation (LDA) for exchange and correlation is the most used 
scheme in Density Functional Theory, but other approximations, which try to correct some 
of its deficiencies (e.g., incapacity to bound extra electrons), may be implemented. One 
of the most popular is a semi-local approach named Generalized Gradient Approximation 



(GGA) from which there are a couple of versions, e.g. Refs. [31] and |32]. Instead of LDA, 
we have used the GGA, in the form recently proposed by Perdew, Burk and Ernzerhof [32 1, 
for calculating the cohesive energies for sodium and aluminum clusters in the SJM as well as 
in the SAPS. In the SAPS we have taken the ionic configuration optimized within the LDA. 
For SJM we are not considering self-compression. 

The results are shown in Fig. [j], in comparison with LDA theory and experimental data. 
We obtain a small decrease of the cohesive energy for both metals. This fact reiterates our 
conclusion, reached with the LDA, that the SJM is in better agreement with the experimental 
data than the SAPS. The decrease of the cohesive energy, due to gradient corrections, is a 
well- known feature of Density Functional Theory ||31|| . 
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Figure 1: Valence electron densities of neutral and charged octamers of sodium and alu- 
minum, obtained with the SAPS and the SJM. In insert, we represent the radial probability 
density of the neutral clusters, r 2 n(r), in arbitrary units. The horizontal line represents 
the background density of the neutral cluster in the SJM (its radius is Rj e \ = r*A^s, with 
r* = 3.70, for sodium, and r* = 1.92 for aluminum. The vertical arrows on the horizontal 
axis denote the position of ionic shells in the SAPS: in each picture the left arrow marks 
the position of the ionic shell for the neutral cluster, whereas the right one stands for the 
charged case. Each cluster has a single ionic shell, with D 4 ^ symmetry. 

Figure 2: Cohesive energies of sodium and aluminum clusters obtained in the SAPS model 
and the SJM. The SJM is shown in two versions: one, denoted by SJM(rf), with the 
background density fixed at the bulk value, and another denoted by SJM(r*) including self- 
compression, where the background density of the cluster is allowed to change. The solid 
lines are obtained with the liquid drop formula, with (below) and without (above) the self- 
compression effect. (We take the Kohn-Sham energy of the atom instead of the liquid drop 



value.) Ou results are compared with Car-Parrinello calculations [18, 20 



Figure 3: Binding energies per atom for sodium clusters in the SJM and SAPS models. The 
lines marked with TFDGW4 represents the energy obtained by replacing the quantal kinetic 
energy by the TFDGW4 functional of the Kohn-Sham density. The numbers denote the 
structure of ionic shells in SAPS, and have the following meaning: 1- one ionic shell; 2- one 
ionic shell with one atom at the center; 3- two ionic shells; and 4- two ionic shells with one 
atom at the center. 



Figure 4: Size evolution of the outer ionic shell radius, R ou t, in the SAPS, and of the 
equilibrium ionic density of SJM, r*, for sodium clusters and aluminum. Both quantities are 
normalized to the corresponding limit N — > oo. The numbers over the SAPS points have 
the same meaning as in Fig. 3. 



Figure 5: Size of Na§ + and Alg + , normalized to the corresponding size of the neutral system, 
as a function of charge q. For the SAPS we represent the ratio between the outer ionic shell 
radius of the charged and the neutral cluster. For the SJM we represent the ratio between 
the equilibrium density parameter, r*, of the charged and the neutral cluster. The inset is a 
blow-up between q = —2 and q = 3 (for Al| the SAPS has no bound solutions for q < 0). 



Figure 6: Elastic stiffness B of Alg + , as a function of charge q, obtained with the SAPS and 
the SJM. 



Figure 7: Cohesive energies of sodium (above) and aluminum (below) clusters obtained 
within the SAPS and the SJM, with the LDA and GGA for exchange and correlation. 
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